################################################################################
# Libraries and Import
################################################################################
rm(list=ls())

library(sandwich)
library(tidyverse)
library(ri2)
library(stargazer)

load(file="data/dat_analysis_JP.RData") 

set.seed(30) # Randomization inference values will change based on seed

# Control number of simulations in randomization procedures
sims = 10000

################################################################################
# Define heterogenous groups
################################################################################
# Define party as ruling party
vignette$dpj <- with(vignette, ifelse((party == 6), 1, 0))
vignette$party <- with(vignette, ifelse((party == 5 | party == 9), 1, 0))

# Redefine education variable to binary college or no college
vignette$edu <- with(vignette, ifelse((edu == 1 | edu == 2), 0, 1))

# Liberal peace variables
vignette$China_trade <- vignette$China_export + vignette$China_import
vignette$inter_pos <- as.numeric(vignette$international_sys==4) + 
                      as.numeric(vignette$international_sys==5)

# Create age and age-squared variables
vignette$age <- as.numeric(as.character(vignette$age_group))
vignette$age2 <- vignette$age^2

# Nationalism
vignette$perception_fp = as.numeric(as.character(vignette$perception_fp))
vignette$international_sys = as.numeric(as.character(vignette$international_sys))
vignette$territory = as.numeric(as.character(vignette$territory))
vignette$Japan_citizen = as.numeric(as.character(vignette$Japan_citizen))
vignette$critique = as.numeric(as.character(vignette$critique))

dv.names = c("perception_fp", "territory", "Japan_citizen", "critique")
dv.names = c("territory", "Japan_citizen", "critique")
factor.obj = princomp(vignette[, dv.names], cor=TRUE)
vignette$nat_pca = as.vector(factor.obj$scores[,1])

################################################################################
# Split into two datasets by treatment group
################################################################################
# Split dataset by treatment group
nat = vignette %>% filter(!is.na(nat)) %>% rename(z = nat)
econ = vignette %>% filter(!is.na(econ)) %>% rename(z = econ)

################################################################################
# Standard interaction tests
################################################################################
# Nationalist treatment
nat_party = lm(D1 ~ z + party + z:party, data = nat)
nat_education = lm(D1 ~ z + edu + z:edu, data = nat)
nat_trade = lm(D1 ~ z + China_trade + z:China_trade, data = nat)
nat_nationalist = lm(D1 ~ z + nat_pca + z:nat_pca, data = nat)

# Economic treatment
econ_party = lm(D1 ~ z + party + z:party, data = econ)
econ_education = lm(D1 ~ z + edu + z:edu, data = econ)
econ_trade = lm(D1 ~ z + China_trade + z:China_trade, data = econ)

################################################################################
# Export tables with standard interaction test
################################################################################
# Nationalism
stargazer(nat_nationalist, 
          out = "tables/TA6. japan_nat_nationalist.tex", 
          se = starprep(nat_nationalist),
          title="Nationalism and support for nationalist interference (Japan)",
          label = "tab:nat_nationalist",
          align=TRUE, 
          dep.var.labels=c("Government Support"),
          keep = c("z", "nat_pca", "z:nat_pca", "Constant"),
          covariate.labels=c("Government interference", 
                             "Nationalist", 
                             "Government interference x nationalist", 
                             "Constant"),
          no.space=TRUE,
          keep.stat = c("n")
          )

# Liberal peace
stargazer(nat_trade, nat_education, econ_trade, econ_education,
          out = "tables/TA13. japan_edu_trade.tex",
          se = starprep(nat_trade, nat_education, econ_trade, econ_education),
          title="Testing the Liberal Peace Logic (Japan)",
          label = "tab:japan_edu_trade",
          align=TRUE,
          dep.var.labels=c("Government Support"),
          column.labels=c("Nationalist", "Nationalist", "Economic", "Economic"),
          keep = c("z", "China_trade", "z:China_trade", "edu", "z:edu", "Constant"),
          covariate.labels=c("Government interference",
                             "Trade with China",
                             "Government interference x trade",
                             "College education",
                             "Government interference x college",
                             "Constant"),
          no.space=TRUE,
          keep.stat = c("n"))

################################################################################
# Individual CATEs method: Government action support by Party
################################################################################
# Nationalist: Calculate difference in CATEs by LDP support
nat_party_1 = lm(D1 ~ z , data = nat, subset = party == 1)
nat_party_2 = lm(D1 ~ z, data = nat, subset = party == 0)
dic <- abs(nat_party_1$coefficients[[2]] - nat_party_2$coefficients[[2]])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- nat
dic_null <- vector(mode = "integer", length = sims)

# 1. # Form full schedule of potential outcomes assuming constant effects
null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 

# 2. Assign subjects randomly to treatment or control
  null$z_null = complete_ra(nrow(null))

# 3. Calculate difference between estimated CATEs
  nat_party_null_1 = lm(D1_null ~ z_null, data = null, subset = party == 1)
  nat_party_null_2 = lm(D1_null ~ z_null, data = null, subset = party == 0)
  dic_null[[i]] <- abs(nat_party_null_1$coefficients[[2]] - nat_party_null_2$coefficients[[2]])
}

# 4. Calculate two sided p-value (0.06)
p_nat_party <- sum(abs(dic_null) >= dic)/length(dic_null)

############# Economic: Calculate difference in CATEs by LDP support ###########
econ_party_1 = lm(D1 ~ z, data = econ, subset = party == 1)
econ_party_2 = lm(D1 ~ z, data = econ, subset = party == 0)
dic <- abs(econ_party_1$coefficients[2] - econ_party_2$coefficients[2])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- econ
dic_null <- vector(mode = "integer", length = sims)

null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 

  null$z_null = complete_ra(nrow(null))
  
  nat_party_null_1 = lm(D1_null ~ z_null, data = null, subset = party == 1)
  nat_party_null_2 = lm(D1_null ~ z_null, data = null, subset = party == 0)
  
  dic_null[[i]] <- abs(nat_party_null_1$coefficients[2] - nat_party_null_2$coefficients[2])
}

p_econ_party <- sum(abs(dic_null) >= dic)/length(dic_null)

################################################################################
# Government action support by education level
################################################################################
# Nationalist: Calculate difference in CATEs by education (college or no)
nat_edu_1 = lm(D1 ~ z, data = nat, subset = edu == 1)
nat_edu_2 = lm(D1 ~ z, data = nat, subset = edu == 0)
dic <- abs(nat_edu_1$coefficients[2] - nat_edu_2$coefficients[2])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- nat
dic_null <- vector(mode = "integer", length = sims)

null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 

  null$z_null = complete_ra(nrow(null))
  
  nat_edu_null_1 = lm(D1_null ~ z_null, data = null, subset = edu == 1)
  nat_edu_null_2 = lm(D1_null ~ z_null, data = null, subset = edu == 0)
  
  dic_null[[i]] <- abs(nat_edu_null_1$coefficients[2] - nat_edu_null_2$coefficients[2])
}

p_nat_edu <- sum(abs(dic_null) >= dic)/length(dic_null)

############# Economic: Calculate difference in CATEs by education #############
econ_edu_1 = lm(D1 ~ z, data = econ, subset = edu == 1)
econ_edu_2 = lm(D1 ~ z, data = econ, subset = edu == 0)
dic <- abs(econ_edu_1$coefficients[2] - econ_edu_2$coefficients[2])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- econ
dic_null <- vector(mode = "integer", length = sims)

null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 

  null$z_null = complete_ra(nrow(null))
  
  nat_edu_null_1 = lm(D1_null ~ z_null, data = null, subset = edu == 1)
  nat_edu_null_2 = lm(D1_null ~ z_null, data = null, subset = edu == 0)
  
  dic_null[[i]] <- abs(nat_edu_null_1$coefficients[2] - nat_edu_null_2$coefficients[2])
}

p_econ_edu <- sum(abs(dic_null) >= dic)/length(dic_null)

################################################################################
# Government action support by trade relationship
################################################################################
# Nationalist: Calculate difference in CATEs by education (college or no)
nat_trade_1 = lm(D1 ~ z, data = nat, subset = China_trade == 1)
nat_trade_2 = lm(D1 ~ z, data = nat, subset = China_trade == 0)
dic <- abs(nat_trade_1$coefficients[2] - nat_trade_2$coefficients[2])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- nat
dic_null <- vector(mode = "integer", length = sims)

null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 

  null$z_null = complete_ra(nrow(null))
  
  nat_trade_null_1 = lm(D1_null ~ z_null, data = null, subset = China_trade == 1)
  nat_trade_null_2 = lm(D1_null ~ z_null, data = null, subset = China_trade == 0)
  
  dic_null[[i]] <- abs(nat_trade_null_1$coefficients[2] - nat_trade_null_2$coefficients[2])
}

p_nat_trade <- sum(abs(dic_null) >= dic)/length(dic_null)

############# Economic: Calculate difference in CATEs by education #############
econ_trade_1 = lm(D1 ~ z, data = econ, subset = China_trade == 1)
econ_trade_2 = lm(D1 ~ z, data = econ, subset = China_trade == 0)
dic <- abs(econ_trade_1$coefficients[2] - econ_trade_2$coefficients[2])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- econ
dic_null <- vector(mode = "integer", length = sims)

null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 

  null$z_null = complete_ra(nrow(null))
  
  nat_trade_null_1 = lm(D1_null ~ z_null, data = null, subset = China_trade == 1)
  nat_trade_null_2 = lm(D1_null ~ z_null, data = null, subset = China_trade == 0)
  
  dic_null[[i]] <- abs(nat_trade_null_1$coefficients[2] - nat_trade_null_2$coefficients[2])
}

p_econ_trade <- sum(abs(dic_null) >= dic)/length(dic_null)

################################################################################
# Government action support by pre-existing nationalism
################################################################################
# Nationalist: Calculate difference in CATEs by nationalism
nat_nat_1 = lm(D1 ~ z, data = nat, subset = Japan_citizen == 1)
nat_nat_2 = lm(D1 ~ z, data = nat, subset = Japan_citizen == 0)

nat_nat_1 = lm(D1 ~ z, data = nat, subset = nat_pca > mean(nat_pca))
nat_nat_2 = lm(D1 ~ z, data = nat, subset = nat_pca < mean(nat_pca))
dic <- abs(nat_nat_1$coefficients[2] - nat_nat_2$coefficients[2])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- nat
dic_null <- vector(mode = "integer", length = sims)

null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 

  null$z_null = complete_ra(nrow(null))
  
  #nat_nat_null_1 = lm(D1_null ~ z_null, data = null, subset = Japan_citizen == 1)
  #nat_nat_null_2 = lm(D1_null ~ z_null, data = null, subset = Japan_citizen == 0)
  
  nat_nat_null_1 = lm(D1_null ~ z_null, data = null, subset = nat_pca < mean(nat_pca))
  nat_nat_null_2 = lm(D1_null ~ z_null, data = null, subset = nat_pca < mean(nat_pca))
  
  dic_null[[i]] <- abs(nat_nat_null_1$coefficients[2] - nat_nat_null_2$coefficients[2])
}

p_nat_nat <- sum(abs(dic_null) >= dic)/length(dic_null)

################################################################################
# Look at DPJ supporters in Japan: exploratory only
################################################################################
########### Nationalist: Calculate difference in CATEs by DPJ support ##########
nat_party_1 = lm(D1 ~ z, data = nat, subset = dpj == 1)
nat_party_2 = lm(D1 ~ z, data = nat, subset = dpj == 0)
dic <- abs(nat_party_1$coefficients[[2]] - nat_party_2$coefficients[[2]])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- nat
dic_null <- vector(mode = "integer", length = sims)

# 1. # Form full schedule of potential outcomes assuming constant effects
null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 
  
  # 2. Assign subjects randomly to treatment or control
  null$z_null = complete_ra(nrow(null))
  
  # 3. Calculate difference between estimated CATEs
  nat_party_null_1 = lm(D1_null ~ z_null, data = null, subset = party == 1)
  nat_party_null_2 = lm(D1_null ~ z_null, data = null, subset = party == 0)
  dic_null[[i]] <- abs(nat_party_null_1$coefficients[[2]] - nat_party_null_2$coefficients[[2]])
}

# 4. Calculate two sided p-value (0.06, cannot reject null of constant effects,
# especially once Bonferroni correction is taken into account)
p_nat_dpj <- sum(abs(dic_null) >= dic)/length(dic_null)

############# Economic: Calculate difference in CATEs by DPJ support ###########
econ_party_1 = lm(D1 ~ z, data = econ, subset = dpj == 1)
econ_party_2 = lm(D1 ~ z, data = econ, subset = dpj == 0)
dic <- abs(econ_party_1$coefficients[2] - econ_party_2$coefficients[2])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- econ
dic_null <- vector(mode = "integer", length = sims)

null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 
  null$z_null = complete_ra(nrow(null))
  
  nat_party_null_1 = lm(D1_null ~ z_null, data = null, subset = party == 1)
  nat_party_null_2 = lm(D1_null ~ z_null, data = null, subset = party == 0)
  
  dic_null[[i]] <- abs(nat_party_null_1$coefficients[2] - nat_party_null_2$coefficients[2])
}

p_econ_dpj <- sum(abs(dic_null) >= dic)/length(dic_null)

################################################################################
# Calculate adjusted p values
################################################################################
options(scipen=999)

Treatment = c("Nationalist", "Nationalist", "Nationalist", "Nationalist", 
              "Economic", "Economic", "Economic")

Covariate = c("Party", "Education", "Trade", "Nationalism",
              "Party", "Education", "Trade")

Unadjusted = c(p_nat_party, p_nat_edu, p_nat_trade, p_nat_nat, p_econ_party, 
               p_econ_edu, p_econ_trade)

Unadjusted = round(Unadjusted, 4)

Bonferroni = round(p.adjust(Unadjusted, method = "bonferroni"), 4)

BH = round(p.adjust(Unadjusted, method = "BH"), 4)

Holm = round(p.adjust(Unadjusted, method = "holm"), 4)

corrections = cbind(Treatment, Covariate, Unadjusted, BH, Holm, Bonferroni)

# Output table
stargazer(corrections, 
          header=FALSE, 
          summary = FALSE,
          column.sep.width = "2",
          title = 'Multiple comparisons corrections Japan',
          label = "corrections_japan",
          out = "tables/TA9. corrections_japan.tex")

################################################################################
# Liberal peace theory exploratory analysis
################################################################################
nat_model_trade = lm(D1 ~ z + China_trade + z:China_trade, data = nat)
nat_model_trade_cov = lm(D1 ~ z + prov + age + age2 + male + edu + China_trade + z:China_trade, data = nat)

nat_model_influence = lm(D1 ~ z + China_influence + z:China_influence, data = nat)
nat_model_influence_cov = lm(D1 ~ z + prov + age + age2 + male + edu + China_influence + z:China_influence, data = nat)

nat_model_inter = lm(D1 ~ z + inter_pos + z:inter_pos, data = nat)
nat_model_inter_cov = lm(D1 ~ z + prov + age + age2 + male + edu + inter_pos + z:inter_pos, data = nat)

# Output table using stargazer
stargazer(nat_model_trade_cov, nat_model_influence_cov, nat_model_inter,
          out = "tables/TA15. japan_liberal_peace.tex", 
          title="Testing the Liberal Peace Logic (Japan)",
          label = "libpeacejapan",
          align=TRUE, 
          dep.var.labels=c("Government Support", "Government Support", "Government Support"),
          keep = c("z", "China_trade", "z:China_trade", "China_influence", 
                   "z:China_influence", "inter_pos", "z:inter_pos", "Constant"),
          covariate.labels=c("Government interference", 
                             "Trade with China", 
                             "Government interference x Trade with China", 
                             "Chinese products",
                             "Government interference x Chinese products",
                             "Internationalist",
                             "Government interference x Internationalist",
                             "Constant"),
          no.space=TRUE,
          keep.stat = c("n","rsq", "adj.rsq"),
          add.lines = list(c("Covariate Adjustment", "Yes", "Yes", "Yes")))